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Abstract-  A  new  adaptive  signal-preserving  technique  for 
noise  suppression  in  functional  magnetic  resonance  imaging 
(fMRI)  data  is  proposed  based  on  spectrum  subtraction.  The 
proposed  technique  estimates  a  model  for  the  power  spectrum  of 
random  noise  from  the  acquired  data.  This  model  is  used  to 
estimate  a  noise-suppressed  power  spectrum  for  any  given  pixel 
time  course  by  simple  subtraction  of  power  spectra.  The  new 
technique  is  tested  using  computer  simulations  and  real  data  for 
event-related  fMRI  experiments.  The  results  show  the  potential 
of  the  new  technique  in  suppressing  noise  while  preserving  the 
other  deterministic  components.  Moreover,  further  analysis 
using  principal  component  analysis  (PCA)  and  independent 
component  analysis  (ICA)  shows  a  significant  improvement  in 
both  convergence  and  clarity  of  results  when  the  new  technique 
is  used.  This  suggests  the  value  of  the  new  technique  as  a  useful 
preprocessing  step  for  this  type  of  signals. 
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I.  Introduction 

Functional  Magnetic  Resonance  Imaging  (fMRI)  provides 
a  valuable  noninvasive  tool  for  investigating  brain  function.  It 
localizes  brain  activity  during  mental  or  physical  activity  by 
detecting  the  corresponding  increase  in  average  cerebral 
blood  oxygenation  or  cerebral  blood  flow  [1].  To  observe 
these  hemodynamic  changes,  rapid  acquisition  of  a  series  of 
brain  images  is  performed.  The  sequence  of  images  is 
analyzed  to  detect  such  changes  and  the  result  is  expressed  in 
the  form  of  a  map  of  the  activated  regions  in  the  brain. 

Classically,  most  of  fMRI  studies  were  conducted  using 
the  so-called  block  design  approach,  whereby  two  sets  of  data 
are  acquired.  First,  a  number  of  frames  are  acquired  while  the 
subject  is  at  rest  or  under  some  baseline  condition,  then 
another  set  is  acquired  during  the  stimulus  [1].  This  pattern  is 
repeated  for  a  number  of  cycles  in  order  to  improve  SNR, 
which  would  otherwise  be  quite  low.  Recent  advances  in  both 
data  acquisition  and  analysis  have  improved  the  temporal 
resolution  of  fMRI  and  made  it  possible  to  observe  transient 
hemodynamic  changes  with  reasonable  accuracy.  A  good 
example  for  that  is  a  new  experimental  design,  similar  to  that 
of  evoked-response  potential  (ERP)  protocol,  called  single 
trial  or  event-related  fMRI.  In  this  new  design,  the  subject 
receives  a  short  stimulus  or  performs  a  single  instance  task 
while  the  resultant  transient  response  is  measured  [2].  ER- 
fMRI  offers  many  advantages  over  block  design  that  include 
versatility,  investigation  of  trial-to-trial  variations,  extraction 
of  epoch-dependent  information  and  direct  adaptation  of  the 
methods  used  for  ERP  to  fMRI  [2].  The  main  drawback  of 
ER-fMRI  is  the  degradation  in  signal  to  noise  ratio  (SNR)  due 
to  the  transient  nature  of  the  response.  As  a  result,  such 
studies  now  include  epoch  averaging.  Nevertheless,  this 
comes  at  the  expense  of  suppressing  the  information  about 
intra-subject  variations  related  to  psychophysiological 
function  with  each  execution  of  the  task.  Therefore,  a 


processing  method  that  can  be  used  to  suppress  noise  in  the 
acquired  data  would  be  very  useful  to  reduce  the  experiment 
duration  and  preserve  the  information  within  the  acquired 
data. 

Several  methods  of  data  analysis  have  been  used  to 
process  the  fMRI  raw  data.  The  ultimate  goal  of  such  analysis 
is  to  try  to  separate  signal  components  due  to  true  activation, 
physiological  fluctuations  and  random  noise.  The  latter  two 
components  are  considered  as  nuisance  and  must  be  removed 
for  correct  results.  Among  the  most  powerful  techniques  that 
can  be  used  to  separate  signal  components  are  those  based  on 
blind  source  separation  such  as  principal  component  analysis 
(PCA)  and  independent  component  analysis  (ICA).  These 
techniques  decompose  the  signal  sources  into  either 
orthogonal  components  (PCA)  or  more  generally  independent 
components  (ICA).  According  to  the  assumptions  of  both 
techniques,  the  number  of  independent  signal  components 
must  be  less  than  or  equal  to  the  number  of  signals  to  be 
analyzed.  Otherwise,  the  separation  of  components  yields 
incorrect  results  or  even  may  not  converge  at  all  as  in  ICA. 
Unfortunately,  this  condition  is  not  satisfied  in  fMRI  data 
sets.  Given  the  general  assumption  of  uncorrelated  noise,  the 
number  of  components  of  random  noise  alone  is  equal  to  the 
number  of  signals.  The  total  number  of  components  has  to 
add  the  number  of  components  due  to  physiological 
fluctuations  as  well  as  the  activation  components.  As  a  result, 
the  use  of  PCA  and  ICA  based  techniques  may  not  yield 
useful  results  in  this  case.  Therefore,  a  technique  that 
suppresses  random  noise  or  removes  some  of  its  components 
would  be  rather  useful  for  making  the  use  of  PCA  and  ICA 
more  robust  for  clinical  practice. 

In  this  work,  we  study  the  problem  of  reducing  the 
random  noise  while  preserving  the  other  deterministic 
components  in  fMRI  signals.  A  new  adaptive  technique  is 
proposed  based  on  spectrum  subtraction.  The  theoretical 
analysis  of  the  new  technique  and  the  implementation  details 
are  presented.  The  new  technique  is  tested  using  computer 
simulations  as  well  as  real  data  and  the  performance  is 
analyzed.  Finally,  the  value  of  the  proposed  method  as  a 
preprocessing  stage  for  PCA/ICA  techniques  is  demonstrated. 

II.  Theory 

Generally  speaking,  the  fMRI  temporal  signal  can  be 
modeled  as  the  summation  of  the  true  activation  signal,  a 
physiological  baseline  fluctuation  component,  and  a  random 
noise  component.  The  physiological  baseline  fluctuation 
component  can  be  considered  as  a  deterministic  yet  unknown 
signal.  Therefore,  we  will  consider  a  model  that  is  composed 
of  the  sum  of  one  deterministic  component  d(.)  incorporating 
both  the  true  signal  and  the  physiological  noise  and  an 
uncorrelated  stochastic  component  n(.).  That  is, 
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SCO  =  d(0  +  n(t).  (1) 

Since  these  two  component  are  assumed  independent,  the 
corresponding  power  spectrums  are  related  by, 

PssO)=  Pdd(^)+  Pnn(^),  (2) 

where  cross  terms  vanish  because  the  two  components  are 
assumed  uncorrelated.  Hence,  an  estimate  of  the  power 
spectrum  of  the  deterministic  component  takes  the  form, 

PddO)=  PssO)  -  PnnO).  (3) 

That  is,  the  signal  power  spectrum  is  obtained  by  spectrum 
subtraction  of  the  noisy  signal  and  noise  power  spectra.  In 
order  to  compute  the  deterministic  signal  component  from  its 
power  spectrum,  the  magnitude  of  the  Fourier  transform  can 
be  obtained  as  the  square  root  of  the  power  spectrum.  The 
problem  now  becomes  that  of  reconstructing  the  signal  using 
magnitude  only  information  about  its  Fourier  transform. 
Several  techniques  can  be  used  to  do  that.  The  one  used  for 
this  work  relies  on  an  estimate  obtained  from  the  phase  of  the 
Fourier  transform  of  the  original  signal.  Hence,  the  Fourier 
transform  of  the  processed  signal  can  be  expressed  as, 

De(t)=  Pdd(<y)1/2 .  Exp(j  Phase(S(<y))).  (4) 

The  enhanced  deterministic  signal  is  just  the  inverse  Fourier 
transform  of  this  expression. 

III.  Methods 

A.  Adaptive  Noise  Model  Estimation 

In  fMRI,  the  acquired  data  set  usually  contains  large  areas 
of  background  and  areas  without  activation.  The  time  courses 
of  pixels  within  such  areas  can  be  used  to  estimate  a  suitable 
model  for  the  noise  power  spectrum.  In  our  implementation, 
only  background  pixels  (defined  by  simple  intensity  masking 
of  the  images)  were  used.  Two  methods  can  be  used  to 
compute  an  estimate  of  the  power  spectrum  of  noise  using 
parametric  and  nonparametric  approaches. 

In  the  parametric  approach,  a  noise  model  is  assumed  and 
the  model  parameters  are  estimated  from  the  data.  The  power 
spectrum  is  subsequently  calculated  from  the  model.  In  the 
case  of  fMRI  data,  such  model  can  be  assumed  as  a  zero- 
mean  white  Gaussian  noise.  Consequently,  the  power 
spectrum  can  be  simply  obtained  as  a  flat  curve  with 
magnitude  equal  to  an  estimate  of  the  variance  of  the 
background  areas.  On  the  other  hand,  the  nonparametric 
approach  does  not  assume  a  particular  model  for  noise.  The 
averaged  periodogram  estimate  for  the  noise  power  spectrum 
is  obtained  directly  from  the  pixel  time  courses  of 
background  areas.  Since  the  number  of  pixels  in  such  areas  is 
expected  to  be  large,  the  variance  of  such  estimate  is  expected 
to  be  very  small. 

The  difference  between  the  two  estimation  approaches  is 
that  the  parametric  technique  models  the  Gaussian  random 
noise  component  of  the  original  signal,  while  the 
nonparametric  technique  may  also  include  other  components 
such  as  global  baseline  variations.  The  selection  of  which  to 
use  depends  on  the  type  of  subsequent  analysis.  In  this  work, 
we  present  the  results  from  the  parametric  approach  to  make 
the  analysis  of  the  results  consistent  in  comparison  and  to 
keep  the  baseline  variations  in  the  processed  data  to  assess  the 
performance  of  PCA  and  ICA  in  estimating  such  components. 


B.  Signal  Power  Spectrum  Estimation 

Since  the  proposed  technique  is  applied  to  a  single  time 
course  at  a  time,  the  periodogram  estimate  of  signal  power 
spectrum  is  expected  to  have  a  rather  large  variance  [3].  As  a 
result,  the  subtraction  of  power  spectra  in  (3)  may  contain 
negative  values  in  practical  implementations.  This  causes  a 
problem  in  trying  to  compute  the  square  root  to  recover  the 
processed  signal.  A  simple  approach  to  overcome  this 
problem  is  to  replace  all  negative  values  in  the  subtraction 
results  by  zero  [3]. 

C.  Statistical  Noise  Removal 

Given  the  nature  of  the  original  signal,  we  observe  that 
the  variance  in  the  power  spectrum  estimate  may  only  result 
from  the  random  component.  Since  the  expected  value  of  the 
noise  variation  is  known  from  the  derived  model  and  given 
the  statistical  characteristics  of  the  periodogram  estimate,  we 
can  express  the  noise  at  each  of  the  power  spectrum 
frequency  bins  as  a  Gaussian  random  variable  with  mean  and 
variance  both  equal  to  the  noise  model  [3].  As  a  result,  the 
subtraction  in  (3)  would  effectively  remove  only  a  part  of  the 
noise  power  spectrum.  In  other  words,  the  upper  half  of  the 
Gaussian  distribution  would  still  remain  in  the  processed 
signal. 

To  solve  this  problem,  a  slight  modification  to  the 
technique  is  added  to  allow  direct  control  over  the  extent  of 
noise  removed.  The  modified  equation  takes  the  form, 

Pdd(^)=  Pss(^)  -  a  .  Pnn(^).  (5) 

Here,  the  factor  a  is  added  to  control  the  confidence  of  noise 
removal.  This  problem  can  be  expressed  in  the  form  of  a 
statistical  z-test  where  the  threshold  a  controls  the  p-value  of 
the  test.  That  is,  the  larger  the  value  of  a,  the  less  the 
probability  that  the  output  power  spectrum  contains  a  noise 
component.  On  the  other  hand,  increasing  this  value  would 
increase  the  likelihood  that  some  parts  of  the  signal  may  also 
be  removed.  Therefore,  the  selection  of  the  value  of  a  is 
useful  to  fine-tune  the  results  of  the  new  technique.  Several 
optimization  criteria  can  be  used  to  select  the  value  of  this 
parameter.  An  example  of  these  is  the  use  of  entropy  based 
objective  function  optimized  over  the  autocorrelation  function 
of  the  difference  between  the  original  and  processed  signals 
for  different  a  values.  This  favors  the  values  of  a  that  give  an 
autocorrelation  function  with  narrow  extent  around  zero  and 
of  minimal  side  peaks.  This  tends  to  preserve  the  components 
of  the  true  signal,  which  give  rise  to  periodic  peaks  in  the 
autocorrelation  function.  In  this  work,  we  used  a  fixed  value 
of  this  parameter  that  is  equal  to  1  to  make  it  easier  to 
compare  the  results  and  assess  the  improvement  after  using 
this  technique  as  a  preprocessing  stage. 

D.  Analysis  of  Results  Using  PCA  and  ICA 

To  show  the  improvement  in  using  PCA  and  ICA  on  the 
processed  signal,  the  new  technique  is  applied  to  process  all 
pixel  time  courses  in  the  acquired  data  set  independently  and 
then  the  processed  data  set  is  used  for  subsequent  PCA/ICA. 
The  PCA/ICA  techniques  were  performed  using  a  Matlab 
(Math  Works,  Inc.)  program  [4].  The  goal  of  this  analysis  is 


to  assess  the  performance  of  the  new  technique  in  enhancing 
the  results  of  PCA  and  ICA  and  stabilizing  the  convergence 
characteristics  of  the  ICA.  Moreover,  the  difference  signals 
between  the  original  and  filtered  data  sets  were  also  analyzed 
using  these  techniques.  This  helps  verify  the  absence  of  signal 
components  within  this  discarded  part  of  the  original  signal. 

IV.  Experimental  verification 

The  proposed  technique  was  verified  using  computer 
simulations  as  well  as  actual  data  from  a  human  volunteer. 
The  computer  simulations  were  performed  in  a  similar 
fashion  to  [2]  whereby  a  computer  generated  event-related 
fMRI  activation  signal  was  added  to  an  actual  baseline  data 
set.  The  baseline  data  were  collected  on  a  healthy  human 
volunteer  using  an  EPI  sequence  (TE/TR=25/500  ms,  FOV= 
20cm  x  20cm,  slice  thickness=5mm,  images  640)  on  a 
Siemens  Magnetom  Vision  1.5T  clinical  scanner.  The  number 
of  epochs  was  8  and  the  length  of  each  epoch  was  64.  The 
generated  activation  was  designed  to  include  inter-epoch 
variations  in  both  the  magnitude  and  width  of  the  activation 
signal  in  order  to  test  the  performance  of  the  new  technique  in 
preserving  such  variations.  The  overall  standard  deviation  of 
the  generated  activation  was  varied  to  test  the  performance 
under  different  values  for  the  signal-to-noise  ratio  (SNR)  [2]. 

The  actual  data  were  obtained  from  an  event-related  fMRI 
study  performed  on  a  normal  human  volunteer  using  a 
Siemens  1.5T  Magnetom  Vision  clinical  scanner  [2].  In  this 
study,  an  oblique  slice  through  the  motor  and  the  visual 
cortices  was  imaged  using  a  T2  -weighted  EPI  sequence 
(TE/TR=60/300ms,  flip  angle=55°,  FOV=22cm  x  22cm,  slice 
theckness=5mm).  The  subject  performed  rapid  finger 
movement  cued  by  flashing  LED  goggles.  The  study  consists 
of  31  epochs  with  64  images  per  epoch.  Temporal  data  from 
only  8  epochs  of  pixels  in  both  the  motor  and  visual  cortices 
were  processed  using  the  new  method  and  compared  to  the 
average  of  all  epochs.  The  PCA  and  ICA  techniques  were 
applied  to  decompose  the  signal  into  its  basic  components. 
Both  techniques  were  used  to  process  time  course  signals 
before  and  after  the  new  technique  is  applied  on  pixels  within 
a  window  selected  by  the  user.  Moreover,  the  difference 
signals  were  also  analyzed  from  the  same  window. 

V.  Results  and  Discussion 

The  results  of  applying  the  new  technique  to  process 
computer  simulated  fMRI  data  are  shown  in  Fig.  1  for  SNR 
values  of  1.0  and  0.5.  As  can  be  observed,  the  noise  in  the 
original  data  was  suppressed  significantly  in  the  output  and 
the  difference  signal  appears  free  of  signal  components.  In 
Fig.  2,  the  proposed  technique  is  used  to  process  real  fMRI 
data  from  two  activated  pixels.  The  results  look  dramatically 
improved  compared  to  the  original.  In  fact,  they  appear  even 
less  noisy  than  the  result  from  averaging  while  maintaining 
the  signal  structure  unaltered.  We  also  notice  that  the  baseline 
variations  also  remained  unaltered  as  a  result  of  the  use  of  the 
parametric  approach  for  our  noise  model.  Again,  the 
difference  signals  appear  to  have  no  signal  components.  In 
Figs.  3  and  4,  the  results  of  applying  the  PCA  and  ICA 
techniques  before  and  after  processing  with  the  new 


technique  are  presented.  The  figures  show  the  first  four 
components  of  each.  The  results  appear  very  noise  before 
applying  the  new  technique  in  both  methods.  The  results  after 
it  was  used  appear  significantly  clearer.  We  notice  that  the 
baseline  variation  component  now  appears  in  the  ICA  results, 
which  was  not  present  in  the  analysis  before  denoising. 
Moreover,  severe  instability  was  observed  when  using  the 
ICA  iteration  on  noise  data.  In  fact,  the  analysis  of  the 
difference  signal  using  this  method  never  reached 
convergence  in  any  of  our  experiments.  This  supports  our 
hypothesis  of  the  need  to  remove  noise  components  to  make 
the  number  of  components  less  than  the  number  of  signals. 
Also,  this  shows  that  the  proposed  method  was  indeed 
successful  to  remove  such  components  without  affecting  the 
true  signal  form. 

The  results  of  using  the  new  technique  suggest  that  it  has 
the  effect  of  robust  removal  of  random  noise  while  preserving 
deterministic  signal  components.  Because  it  relies  on 
subtracting  the  noise  component,  it  does  not  affect 
independence  of  data  points  within  the  time  course.  The 
technique  also  has  a  very  small  computational  complexity.  It 
also  has  the  advantage  of  being  adaptive  with  very  few 
assumptions  about  the  noise  model  and  no  assumptions  about 
signal.  Finally,  it  provides  statistical  control  over  noise 
removal  using  a  single  parameter.  The  results  indicate  that  the 
new  method  enables  robust  use  of  PCA  and  ICA.  The 
limitations  of  the  new  technique  include  the  fact  that  noise  in 
the  phase  component  is  not  removed  and  that  spectrum 
clipping  may  occur  due  to  large  variance  in  power  spectrum 
estimate  from  periodogram  without  averaging. 

VI.  Conclusions 

A  new  signal  denoising  technique  was  proposed  for  fMRI 
signals.  The  method  is  adaptive  and  simple  to  implement 
while  offering  a  substantial  improvement  of  signal  to  noise 
ratio.  The  new  technique  was  demonstrated  for  computer 
simulated  and  real  data  and  also  shown  to  improve  the 
performance  of  PCA  and  ICA  in  analyzing  fMRI  data  sets. 
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Figure  1:  Results  from  simulations  at  different  SNR  values. 


Figure  2:  Results  from  actual  data  for  two  pixel  time  courses. 
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Figure  3:  PCA/ICA  results  before  processing. 
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Figure  4:  PCA/ICA  results  after  processing. 


